######## INFO ########

# PROJECT
# Replication files for: Economic Impact of the South China Sea Dispute in China-Philippines
#                        Relations from 2012 to 2016: A DM-LFM analysis  

# Authors: Dianyi Yang
# R Script
# Purpose: This script runs the DM-LFM and SCM analysis for the exports to China, produces basic  
#          tables and figures, and store the results for more comprehensive analysis in the next scripts.
# Created: 28 Oct 2023
# Updated: 29 Nov 2023
# Inputs: input/export_to_China.csv
#         custom_functions/RMSPE.r
# Outputs: Table 1 (summary statistics)
#          stored/dmlfm_placebo.rdata
#          stored/dmlfm.rdata
#          stored/leave-one-out.rdata
#          output/sc_trend.pdf
#          output/sc_trend.png
#          sc_RMSPE.pdf
#          sc_RMSPE.png
#          stored/sc.rdata


######## SETUP ########
rm(list = ls()) # clear workspace

need <- c("reshape2","devtools", "synthdid",'tidyverse',
          'bpCausal', "fixest","modelsummary",
          'kableExtra') # list packages needed
have <- need %in% rownames(installed.packages()) # checks packages you have
if(any(!have)) install.packages(need[!have]) # install missing packages
#devtools::install_github("synth-inference/synthdid") #manually install sdid if the codes above fail to do so automatically
invisible(lapply(need, library, character.only=T)) # load needed packages

setwd(gsub('/code','',dirname(rstudioapi::getSourceEditorContext()$path))) #set wd to root directory
list.files()
data <- read.csv("input/export_to_China.csv")
seed=666 #set seed
set.seed(seed)

#logimport
data$logimport <- log(data$C.import)
data$treated <- ifelse(data$Unit == 1 & data$M.unit >= 0, 1,0 ) #data cleaning

##################Table 1 (summary stats)###################
data$Country <- factor(data$Country)
data$Country <- relevel(data$Country, ref="C-Philippines")
datasummary( Country ~ ((`Exports to China` = `C.import`)+(`Log(Exports from China)` = logimport))
             *(Mean + SD + Min + Max),
            data = data,
           # output='latex' #if you want to export to latex
            )

##################time placebo###############################
data <- read.csv("input/export_to_China.csv")
data_placebo <- subset(data, M.unit<0)
data_placebo$logimport <- log(data_placebo$C.import)
data_placebo$treated <- ifelse(data_placebo$Unit == 1 & data_placebo$M.unit >= -24, 1,0 ) #data cleaning for synthdid
set.seed(666)
out = bpCausal(data = data_placebo,
               Yname = "logimport",
               Dname = "treated",
               Xname = c(),
               Zname = c(),
               Aname = c(),
               index = c("Country", "M.unit"),
               re = "both",
               ar1 = TRUE,
               r = 10,
               niter = 15000,
               burn = 5000,
               xlasso = 1,
               zlasso = 1,
               alasso = 1,
               flasso = 1,
               a1 = 0.001,
               a2 = 0.001,
               b1 = 0.001,
               b2 = 0.001,
               c1 = 0.001,
               c2 = 0.001,
               p1 = 0.001,
               p2 = 0.001)  #DM-LFM Warning: a bit time consuming
save(out, data_placebo, file='stored/dmlfm_placebo.rdata')  #save model (you can load it next time)

##################DM-LFM####################################
data <- read.csv("input/export_to_China.csv")
#logimport
data$logimport <- log(data$C.import)
data$treated <- ifelse(data$Unit == 1 & data$M.unit >= 0, 1,0 ) #data cleaning
set.seed(666)
out = bpCausal(data = data,
               Yname = "logimport",
               Dname = "treated",
               Xname = c(),
               Zname = c(),
               Aname = c(),
               index = c("Country", "M.unit"),
               re = "both",
               ar1 = TRUE,
               r = 10,
               niter = 15000,
               burn = 5000,
               xlasso = 1,
               zlasso = 1,
               alasso = 1,
               flasso = 1,
               a1 = 0.001,
               a2 = 0.001,
               b1 = 0.001,
               b2 = 0.001,
               c1 = 0.001,
               c2 = 0.001,
               p1 = 0.001,
               p2 = 0.001)  #DM-LFM Warning: a bit time consuming
save(out, data, file='stored/dmlfm.rdata')  #save model (you can load it next time)

#leave-one-out
set.seed(666)
result.use <- effSummary(out) #extract effect information
plot <- result.use[['est.eff']]
plot$miss_unit <- 'All'
for (x in 2:9){
  data_temp <- subset(data, Unit!=x)
  out_temp = bpCausal(data = data_temp,
                 Yname = "logimport",
                 Dname = "treated",
                 Xname = c(),
                 Zname = c(),
                 Aname = c(),
                 index = c("Country", "M.unit"),
                 re = "both",
                 ar1 = TRUE,
                 r = 10,
                 niter = 15000,
                 burn = 5000,
                 xlasso = 1,
                 zlasso = 1,
                 alasso = 1,
                 flasso = 1,
                 a1 = 0.001,
                 a2 = 0.001,
                 b1 = 0.001,
                 b2 = 0.001,
                 c1 = 0.001,
                 c2 = 0.001,
                 p1 = 0.001,
                 p2 = 0.001)  #DM-LFM Warning: a bit time consuming
  result.use_temp <- effSummary(out_temp) #extract effect information
  plot_temp <- result.use_temp[['est.eff']]
  plot_temp$miss_unit <- paste0('w/o ',gsub('C-','',data$Country[data$Unit==x][1]))
  plot <- rbind(plot, plot_temp)
}

save(plot, file='stored/leave-one-out.rdata')

#Synthetic control
setup = panel.matrices(data, unit = "Country", time = "M.unit", outcome = "logimport", treatment = "treated")
tau.hat.sc = sc_estimate(setup$Y, setup$N0, setup$T0)
source("custom_functions/RMSPE.r")
source("custom_functions/solver.r")

plot(tau.hat.sc,
     treated.name = 'Observed',
     control.name = 'Counterfactual',
     line.width   = 0.75,
     effect.alpha = 0.5,
     onset.alpha  = 1)
ggsave('output/sc_trend.pdf', width = 10, height=5) #Figure 13
ggsave('output/sc_trend.png', width = 10, height=5)

#placebo permutation
placebo_result <- data.frame('Country'=character(), RMSPE.ratio=double())
for (n in 1:9){
  data_temp <- subset(data, select=-treated)
  data_temp$treated <- 0
  data_temp$treated[data_temp$Unit==n & data_temp$M.unit>=0] <- 1
  setup = panel.matrices(data_temp, unit = "Country", time = "M.unit", outcome = "logimport", treatment = "treated")
  tau.sc.temp <- sc_estimate(setup$Y, setup$N0, setup$T0)
  placebo_result <- rbind(placebo_result,
                          data.frame('Country' = unique(data_temp$Country[data_temp$treated==1]),
                                     'RMSPE.ratio' = RMSPE_ratio(tau.sc.temp)))
}

placebo_result <- placebo_result %>% arrange(desc(RMSPE.ratio))

ggplot(placebo_result, aes(x=reorder(Country, RMSPE.ratio), y=RMSPE.ratio))+geom_point(shape=15) +coord_flip()+theme_bw()+
  theme(panel.grid.major.x = element_blank(),panel.grid.minor.x = element_blank())+
  labs(x="Countries", y="Postperiod RMSPE/Preperiod RMSPE")
ggsave('output/sc_RMSPE.pdf', width = 10, height=5) #Figure 15
ggsave('output/sc_RMSPE.png', width = 10, height=5)

#placebo CI and p-value
se.sc <- drop(sqrt(vcov(tau.hat.sc, method='placebo')))
pnorm(tau.hat.sc[[1]]/se.sc, lower.tail=FALSE)#p-value
sc_estimate <- list(estimate = drop(tau.hat.sc[[1]]),
                    se = se.sc,
                    p=drop(pnorm(tau.hat.sc[[1]]/se.sc, lower.tail=FALSE)) )
class(sc_estimate) <- 'sc_estimate'
save(tau.hat.sc, sc_estimate, file='stored/sc.rdata')
